Prototype for integration of descriptive and analysis model


This shell scripts execute different types of analysis codes.
Author

Yuta Nakajima

Published

Jun 8, 2024, 05:22:58

1 Workflow

  graph LR;
      M[OML]-->S[SPARQL];
      S-->API;
      API-->P[Python];
      P-->O[Oml2Owl];
      API-->MATLAB;
      MATLAB-->O;
      API-->MODELICA;
      MODELICA-->O;
      O-->M;

Below chart shows an actual workflow in this prototype.

flowchart LR
    id1(Build) --> id2(Query:SPARQL)
    id2 --> id3[massRollup:R]
    id3 --> id4[deltaV:Python]
    id4 --> id5[fuelMass:R]
    id5 --> id6(owl2oml)

2 First Prototype of Process and Data Flow

This prototype uses a Quarto notebook to orchestrate the different Like openMBEE, Jupyter or Quarto perform the execution of codes. Quarto has a function to execute Python and R together. We can bind the data from Python to R using the reticulate package.

flowchart LR
    id01[(OML)] --> id1
    id1(Build) --> id11[(owl)]
    id11 --> id2(Query:SPARQL)
    id2 --> id3[massRollup:R]
    id3 --> id4[deltaV:Python]
    id4 --> id5[fuelMass:R]
    id5 --> id6(owl2oml)
    id6 --> id01
    id2 -.-> id7((json))
    id7 -.-> id3
    id3 -.-> id8((df_mass))
    id8 -.-> id5
    id2 -.-> id9((orbit))
    id9 -.-> id4
    id4 -.-> id10((dv))
    id10 -.-> id5

In this approach, data is exchanged between different analysis tools like a rugby pass-and-receive.

2.1 load utility files

Code
library(reactable)
library(stringr)

searchDirectory <- function(iteration, pattern, parent_directory){
  for(i in 1:iteration){
      path <- list.files(parent_directory, recursive = TRUE, pattern = pattern, full.names = TRUE)
      if(length(path)){
        return(path)
      }
      parent_directory <- dirname(parent_directory)
  }
  print("file not found")
  return(path)
}

source(searchDirectory(4, "osr_common.R", (getwd())))
source(searchDirectory(4, "massRollupKepler16b.R", (getwd())))
Code
library(reticulate)
use_condaenv("py39")

2.2 Build and Run a Fuseki SPARQL endpoint for oml model query

Code
library(omlhashiR)
## oml_repository <- "../open-source-rover/"
oml_repository <- omlrepo
ret1 <- omlhashiR::oml_refresh()
Code
ret2 <- omlhashiR::oml_stop_Daemon(oml_repository)
Code
ret3 <- omlhashiR::oml_callTask(oml_repository, "build")

2.3 Start Fuseki Server

Code
ret4 <- omlhashiR::oml_callTask(oml_repository, "startFuseki")

2.4 Load in OWL ontologies to Fuseki Server

Code
ret5 <- omlhashiR::oml_callTask(oml_repository, "owlLoad")

2.5 Get Configuration & Mass Properties

Code
library(tansakusuR)
endpoint_url <- "http://localhost:3030/tutorial2/sparql"

configuration_root_iri <- "http://example.com/tutorial2/description/components#orbiter-spacecraft"

2.5.1 SPARQL CODE

Code
query_string <- '
PREFIX mission:     <http://imce.jpl.nasa.gov/foundation/mission#>
PREFIX base:        <http://imce.jpl.nasa.gov/foundation/base#>
PREFIX project:     <http://imce.jpl.nasa.gov/foundation/project#>
PREFIX vim4:        <http://bipm.org/jcgm/vim4#>

SELECT DISTINCT ?c1 ?c1_instancename ?c1_id ?c1_name ?c1_mass ?c2 ?c2_instancename ?c2_id ?c2_name ?c2_mass ?c3 ?c3_instancename ?c3_id ?c3_name ?c4
WHERE {

  VALUES ?c1 { <$configuration_root_iri> }


   OPTIONAL{ 
    ?c1 base:hasIdentifier ?c1_id ;
      base:hasCanonicalName ?c1_name ;
  }
    OPTIONAL {
        ?c1_mass_mag vim4:characterizes ?c1 ;
            vim4:hasDoubleNumber ?c1_mass .
    } 
  OPTIONAL{
    ?c1 base:contains ?c2 ;
    OPTIONAL{
      ?c2 base:hasIdentifier ?c2_id ;
          base:hasCanonicalName ?c2_name .
    }
    OPTIONAL {
        ?c2_mass_mag vim4:characterizes ?c2 ;
            vim4:hasDoubleNumber ?c2_mass .
    }     
    OPTIONAL{
      ?c2 project:isSuppliedBy ?c3 ;
          OPTIONAL{
            ?c3 base:hasIdentifier ?c3_id ;
                base:hasCanonicalName ?c3_name;
                project:isAuthorizedBy ?c4 .
          }          
    }    
  }


  BIND(STRAFTER(STR(?c1), "#") AS ?c1_instancename) .
  BIND(STRAFTER(STR(?c2), "#") AS ?c2_instancename) .
  BIND(STRAFTER(STR(?c3), "#") AS ?c3_instancename) .
 }
ORDER BY ?c2_id

'
query_string <- str_replace(query_string, "\\$configuration_root_iri", configuration_root_iri)

2.5.2 Query

Code
df_query <- send_query(endpoint_url,query_string)

2.6 Binding Variables: toR Parameter Table

  • OML descriptions shall specify this parameter table using IMCE Analysis Vocabularies.
  • OML descriptions shall specify mapping of binding variables for each simulation.
  • Get initial value from oml descriptions by SPARQL Query.

2.6.1 SPARQL CODE

Code
query_param_string <- '
PREFIX base:        <http://imce.jpl.nasa.gov/foundation/base#>
PREFIX mission:        <http://imce.jpl.nasa.gov/foundation/mission#>
PREFIX analysis:        <http://imce.jpl.nasa.gov/foundation/analysis#>
PREFIX sa:        <http://example.com/tutorial2/vocabulary/stateanalysis#>

SELECT DISTINCT ?iri ?statevariable ?value
WHERE {

#  VALUES ?analysisTarget {<http://example.com/tutorial2/description/statedictionary#orbiter-spacecraft.delta-v>}
   VALUES ?analysisTarget {<http://example.com/tutorial2/description/statedictionary#orbiter-spacecraft.mass.wet>}
  VALUES ?componentType { sa:StateVariable }
  VALUES ?relationType { sa:isAffectedBy }

  # Recursive Query to get state variables
  ?analysisTarget sa:isAffectedBy* ?iri ;
#  ?iri a ?componentType
  
  OPTIONAL{
    ?iri sa:hasStateValue ?value;
         }



  BIND(STRAFTER(STR(?iri), "#") AS ?statevariable) .
  BIND(STRAFTER(STR(?targetIri), "#") AS ?target) .
  BIND(STRAFTER(STR(?cpIri), "#") AS ?component) .
 }
ORDER BY ?iri

'

2.6.2 Query

Code
df_query_param <- send_query(endpoint_url,query_param_string)

2.6.3 Binding to R Dataframe

Code
# mapping to parameters to instance
# df_parameters <- data.frame(
#   initOrbit = as.double(df_query_param$value[df_query_param$statevariable=="orbiter-spacecraft.initial-orbit"]),
#   targetOrbit = as.double(df_query_param$value[df_query_param$statevariable=="orbiter-spacecraft.target-orbit"]),
#   m_wet = 0.0,
#   m_dry = 0.0,
#   m_fuel = 0.0,
#   dv = 0.0,
#   I_sp = as.double(df_query_param$value[df_query_param$statevariable=="orbiter-spacecraft.isp"])
# )

mapping <- list(c("dv", "orbiter-spacecraft.delta-v"),
  c("initOrbit", "orbiter-spacecraft.initial-orbit"),
  c("I_sp","orbiter-spacecraft.isp"),
  c("m_dry","orbiter-spacecraft.mass.dry"),     
  c("m_fuel","orbiter-spacecraft.mass.fuel"),
  c("m_wet","orbiter-spacecraft.mass.wet"),
  c("targetOrbit","orbiter-spacecraft.target-orbit")
)

# Convert to dataframe
mapping_df <- mapping %>% 
  map_df(~data.frame(.x[1], .x[2])) 
colnames(mapping_df) <- c("parameter", "statevariable")

df_parameters_before <- left_join(df_query_param, mapping_df, by=c("statevariable"="statevariable")) 

# parameter database for quarto with R
df_parameters <- df_parameters_before %>% 
  mutate(value2 = as.double(replace(value, is.na(value), 0))) %>%
  select(parameter,value2) %>%
  pivot_wider(
#    id_cols = age_cat,
    names_from = parameter,
    values_from = value2
    )
Code
df_table <- df_parameters %>% 
  mutate(id=1) %>%
  pivot_longer(!id, names_to = "parameter", values_to = "value") %>%
  select(-id)

datatable(df_table, options = list(autoWidth = TRUE, pageLength = -1))

2.7 MassRollUp

2.7.1 Tidy Data

Code
df_query$c1_type <- "component"

df_query <- df_query %>%
  mutate(c2_label = paste0(c2_id,": ",c2_instancename)) %>%
  mutate(c1_label = paste0(c1_id,": ",c1_instancename)) 

df_config <- df_query

# just for vis add NA to root
df_vis <- df_config %>%
  add_row(c1=NA,
          c1_instancename=NA,
          c1_id=NA,
          c1_name=NA,
          c1_mass=NA,
          c2=df_config$c1[1],
          c2_instancename=df_config$c1_instancename[1],
          c2_id=df_config$c1_id[1],
          c2_name=df_config$c1_name[1],
          c2_mass=df_config$c1_mass[1],
          c1_type=df_config$c1_type[1],
          c2_label=NA,
          c1_label=df_config$c1_label[1],  
          .before = 1)

plotCollapsibleTreeFromDataframe(df_vis, palette="BluYl", 
                                 parent="c2_label",
                                 child="c1_label",
                                 type="c2_id")

2.7.2 Prepare a Graph Data for mass rollup calculation

Code
library(igraph)

df_g <- df_config %>%
  mutate(parent = c1_instancename) %>%
  mutate(child = c2_instancename) %>%
  mutate(type = c2_instancename) %>%
  # mutate(parent = paste0(p_id, ": ", p_instancename)) %>%
  # mutate(child = paste0(c_id, ": ", c_instancename)) %>%  
  select("parent","child","type")

# df_g <- df_g[-1,]

g <- graph_from_data_frame(df_g, 
                           directed = TRUE, 
                           vertices = NULL)

2.7.3 RollUp Pattern

Code
root <- V(g)[1]

# Depth-first search is an algorithm to traverse a graph. It starts from a root vertex and tries to go quickly as far from as possible.
order <- dfs(g, V(g)[root], order.out = TRUE)$order

df_mel_before <- igraph::as_data_frame(g, what = "vertices") %>%
  arrange(factor(name, levels = names(order)))%>%
  filter(name !="NA")

# ここに、df_mass として massの値を入れたコンフィグを作る必要がある。

df_mass <-  df_config

df_mass <-  df_mass %>%
  mutate(mass = as.numeric(c2_mass)) %>%
  select(c2,
         c2_instancename,
         c2_id,
         c2_name,
         c1_id,
         mass) %>%  
  add_row(c2=df_mass$c1[1],
          c2_instancename=df_mass$c1_instancename[1],
          c2_id=df_mass$c1_id[1],
          c2_name=df_mass$c1_name[1],
          c1_id=NA,
          mass=as.numeric(df_mass$c1_mass[1]),
          .before = 1)%>%
  arrange(c2_id)


colnames(df_mass) <- c("c_iri","c_instancename","c_id","c_name", "p_id", "mass")

df_mel_before <- left_join(df_mel_before, df_mass, by = c("name"="c_instancename"))
Code
namekey="name"
masskey="mass"


df_mass_update <- massRollUp(g, root, df_mel_before, namekey = "name",masskey = "mass") %>%
  select(name, mass)


# リーフ数を取得
df_deg <- data.frame(
  name = names(degree(g)),
  degree = degree(g)-1, # num of components is degree(g)-1
  distance = distances(g)[V(g)[1],]
)

# distanceの値で、componentか、workpackageかを切り分けることができる

df_mel_after <- df_mel_before %>% select(-mass)
df_mel_after <- left_join(df_mel_after, df_mass_update, by=c("name"="name")) 
df_mel_after <- left_join(df_mel_after, df_deg, by=c("name"="name")) 

df_mel_after$componenttype <- ifelse(df_mel_after$distance == 0, "system", 
                            ifelse(df_mel_after$distance == 1, "subsystem", 
                                   ifelse(df_mel_after$distance == 2, "assembly", NA)))
Code
df_table <- df_mel_after %>%
  mutate(totalmass=mass) %>%
  select(c_id, name, totalmass, componenttype) 

df_table$componenttype <- factor(df_table$componenttype)

datatable(df_table, options = list(autoWidth = FALSE, pageLength = -1), filter = list(
  position = 'top', clear = FALSE
))

2.8 Binding Variables : toR

Code
df_parameters$m_dry <- df_mel_after$mass[1]
# df_parameters$m_dry <- 325.47880079

2.9 Binding Variables: toPython

Code
orbit_init=r.df_parameters.initOrbit.values.tolist()[0]
orbit_target=r.df_parameters.targetOrbit.values.tolist()[0]
m_final=r.df_parameters.m_dry.values.tolist()[0]
I_sp=r.df_parameters.I_sp.values.tolist()[0]
# orbit_init=r.initOrbit
# orbit_target=r.targetOrbit
# m_init=1000
# I_sp=r.I_sp

2.10 deltaV ( Using Python)

Code
import os
import sys

print('getcwd:      ', os.getcwd())
getcwd:       /home/runner/work/kepler16b-using-imce-vocabulary/kepler16b-using-imce-vocabulary/src/analysis
Code
sys.path.append(os.getcwd())
sys.path.append(os.getcwd()+'/src/analysis')

# from src.analysis import analysisOrbit
import analysisOrbit
from astropy import units as u

analysisRocket = analysisOrbit.hohmanTransfer(orbit_init=orbit_init, orbit_target=orbit_target)
# analysisRocket = analysisOrbit.hohmanTransfer(orbit_init=r.initOrbit, orbit_target=r.targetOrbit, m_init=r.m_dry, I_sp=r.I_sp)
# analysisRocket = analysisOrbit.hohmanTransfer(orbit_init=400, orbit_target=35786, m_init=5000, I_sp=300)

total_delta_v = analysisRocket.calculate_delta_v()


print(f"Total delta-v: {total_delta_v}")
Total delta-v: 3853.959363102248 m / s

2.11 Binding Variables: toR

Code
df_parameters$dv <- py$total_delta_v /1000.0 # km/s

2.12 fuelMass (Using R)

Code
source(searchDirectory(4, "calcWetMass.R", (getwd())))

# I_sp <- 350
# m_dry <- 325.47880079
# df_parameters$dv <- sum(df_parameters$dv)
# Inverse Problem
df_parameters$m_wet <- calcWetMass(df_parameters$dv, df_parameters$m_dry, df_parameters$I_sp)
df_parameters$m_fuel <- df_parameters$m_wet - df_parameters$m_dry
c(df_parameters$m_wet, df_parameters$m_dry, df_parameters$m_fuel)
[1] 6012.68 1957.00 4055.68

2.13 Finally, Update OML

2.13.1 parameter table updated

Code
df_table <- df_parameters %>% 
  mutate(id=1) %>%
  pivot_longer(!id, names_to = "parameter", values_to = "value") %>%
  select(-id)

datatable(df_table, options = list(autoWidth = TRUE, pageLength = -1))

2.13.2 parameters before and after

Code
df_parameters_after <-  left_join(df_parameters_before, df_table,
                                  by = c("parameter"="parameter"),
                                  suffix=c("_before","_after")) %>%
  mutate(status = case_when(
    as.double(value_before) == value_after ~ "unchanged",
    TRUE ~ "changed"
  ))

2.13.3 sending update query

Code
source(searchDirectory(4, "updateSparqlQuery.R", (getwd())))

df_update <- df_parameters_after %>% 
  filter(status == "changed")

endpoint_url <- "http://localhost:3030/tutorial2-tdb/"

for( i in 1:nrow(df_update) ){
  df <- df_update[i,]
  update_value <- as.character(df$value_after)
  update_iri <- df$iri

  df_ret <- updateSparqlQuery(endpoint_url, update_value, update_iri)
  
  ret <- send_update(endpoint_url = endpoint_url, df_ret$query_string_delete_tdb)
  ret <- send_update(endpoint_url = endpoint_url, df_ret$query_string_insert_tdb)
  
}

# cat(df_ret$query_string_delete_tdb)
# cat(df_ret$query_string_insert_tdb)

2.14 Run Save in OML-Vision

Using OML-Vision GUI, run the save task. Then you can see the oml file is updated.

2.15 Verify Update

This can be verified when you re-load the load task and send a query of query_param_string.

Code
endpoint_url <- "http://localhost:3030/tutorial2/"
df_query <- send_query(endpoint_url, query_param_string)


datatable(df_query, options = list(autoWidth = FALSE, pageLength = -1))

3 Discussion: Proposed Process and Data Flow

Instead of using the passing and recieving approach, alternative way is using the owl fuseki database as a hub to exchange the data. OML descriptions specify the process and mapping of data exchanged between analyssis containers.

To do this, we need to specify the analysis container with

  • input
  • output
  • mapping of input and output with other containers

flowchart LR
    oml[(OML)] -.-> id1
    id1(Build) -.-> owl[(owl)]
    id1 --> id2
    owl -.-> id2(Query:SPARQL)
    id2 --> id3[massRollup:R]
    id3 -.-> owl
    id3 --> id4[deltaV:Python]
    id4 -.-> owl
    id4 --> id5[fuelMass:R]
    id5 -.-> owl
    id5 --> id6(owl2oml)
    owl -.-> id6
    id6 -.-> oml
    classDef codpy fill:#ff0000;
    classDef coder fill:#00ffff;
    class id3,id5 coder;
    class id4 codepy;

4 TODO

df_parameters is not good idea to bind the parameters I think this could be better.

Code
df_parameters2 <- data.frame(
  parameter = c("initOrbit","targetOrbit","I_sp","m_dry","m_fuel","m_wet","dv"),
  value = c(400,35786,350,
          df_parameters$m_dry, 
          df_parameters$m_fuel, 
          df_parameters$m_wet,
          df_parameters$dv),
  type = c("input","input","input","output","input","output","output")
)

# how to access data
df_parameters2$value[df_parameters2$parameter == "m_dry"]
[1] 1957

Here is how to access values. when we add the mapping descriptions to analysis codes, we might automate the process of data exchange.

Code
df = r.df_parameters2
df[df['parameter']=='I_sp'].value.values[0]
350.0

4.1 OPTION: m_fuel ( Using Python codes)

Code
analysisRocket.calculate_init_mass(I_sp, m_final)
print(f"Total delta-v: {total_delta_v}")
Total delta-v: 3853.959363102248 m / s
Code
print(f"Required fuel mass: {analysisRocket.m_fuel}\n final: {analysisRocket.m_final}\n init: {analysisRocket.m_init}")
Required fuel mass: 4055.680381236295
 final: 1957.0
 init: 6012.680381236295

4.2 OPTION: deltaV ( Using R codes)

Code
source(searchDirectory(4, "calcDeltaV.R", (getwd())))

dv <- calcDeltaV(initOrbit=df_parameters$initOrbit, targetOrbit = df_parameters$targetOrbit)
[1] 2.397473 1.456487